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Abstract 

We present up to date investigations of the antiferromagnetic Heisenberg 
icosidodecahedron by means of the Density Matrix Renormalization Group 
method. We compare our results with modern Correlator Product State as 
well as Lanczos calculations. 
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1. Introduction 

Thanks to advanced chemical strategies there exist several chemical re- 
alizations of icosidodecahedral magnetic molecules: Mo72Fe3o[l], W72Fe3o[2] 
(both s = 5/2), Mo 72 Cr 30 [3] (s = 3/2), Mo 72 V 30 [4, 5], and W 72 V 30 [6] (both 
s = 1/2). These molecules are some of the largest magnetic molecules syn- 
thesized to date [7]. Icosidodecahedral magnetic molecules are of special 
interest because they are highly symmetric, frustrated, exist with different 
spin quantum numbers, and are a kind of finite-size version of the Kagome 
lattice antiferromagnet [8] . Fig. 1 shows the structure of the icosidodecahe- 
dron. It is an Archimedian Solid with 12 pentagons and 20 triangles, which 
means that it is geometrically frustrated [9]. 

Experimental investigations of these molecules are in most cases measure- 
ments of the susceptibility as a function of temperature [1, 2, 3, 4, 5, 6] or 
magnetic field [10, 11], or magnetization as a function of the applied magnetic 
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Figure 1: Structure of the icosidodecahedron: the black bullets correspond to the spin 
positions and the lines to interaction paths between them. The right part of the figure 
shows the two-dimensional projection. 

field [1, 11, 12]. These experimental investigations show that the icosidodeca- 
hedral magnetic molecules are antiferromagnetic with a nonmagnetic ground 
state. Other experimental techniques that were applied to these molecules 
are NMR and /zSR [13, 14, 15], INS [16], diffuse (elastic) neutron scattering, 
as well as specific heat measurements [17]. 

These molecules are usually modeled using a simple Heisenberg model 
[1, 2, 3, 4, 6, 16, 17, 18, 19, 20, 21, 22, 23]. Anisotropic terms were con- 
sidered in Refs. [22, 24]. Bond disorder and distortions (i.e., more than 
just one exchange constant in the Heisenberg Hamiltonian) were investigated 
in Refs. [10, 23]. However, since the icosidodecahedral molecules comprise 
N = 30 spins, the numerical exact calculation of T = properties is possible 
only for the s = 1/2 case [8, 25, 26]. Quantum Monte Carlo suffers from the 
negative-sign problem so that small temperatures are not feasible [3, 4, 6]. 
Thermodynamical properties at T > can for s = 1/2 also be calculated 
quasi-exactly using the finite-temperature Lanczos method [21, 23, 27]. For 
s > 1/2, approximations are needed. For s = 3/2 and s = 5/2 systems, the 
classical Heisenberg model was used together with efficient classical Monte 
Carlo algorithms [1, 2, 10, 18, 19, 22, 28]. However, such an approximation 
is inappropriate at very low temperatures. The rotational band approxima- 
tion was used in Refs. [1, 12, 16, 17, 20, 29, 30]. Although being a quantum 
mechanical approximation it misses important features of frustrated systems 
such as magnetization jumps and plateaus or low-lying singlets in the spec- 
trum [25]. Another approximation applied to the icosidodecahedron is spin- 
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wave theory [20, 24]. However, as for the rotational band and the classical 
approximation, it is not clear how accurate this approximation is. 

It has to be emphasized here that although there exist many theoretical 
studies on the icosidodecahedron for s > 1/2, accurate numerical calculation 
for the full Heisenberg model are very rare. The Density Matrix Renor- 
malization Group (DMRG) method allows for treating the full Heisenberg 
Hamiltonian but in a reduced Hilbert space [31, 32]. It relies on a controlled 
truncation of the Hilbert space and allows for the estimation of the accu- 
racy so that it seems to be suited for these systems. In Ref. [33] the DMRG 
method has already been applied to the Heisenberg icosidodecahedron with 
s = 5/2. However, only up to m = 120 density matrix eigenstates were used, 
so that the accuracy of the results is rather limited for such a complicated 
system with a geometry that is not favorable for the DMRG method. 

In this article we apply the DMRG method to the antiferromagnetic 
Heisenberg icosidodecahedron. We focus on the calculation of the lowest 
energies in subspaces of total magnetic quantum number M which allow for 
a calculation of the T = magnetization curve and also gaps which might be 
of importance for spectroscopic methods such as, e.g., Inelastic Neutron Scat- 
tering (INS). These results are compared with very recent variational Monte 
Carlo calculations using correlator product states (CPS) by Neuscamman 
and Chan in Ref. [34]. Finally, for the case s = 1/2 we also calculate the dy- 
namical correlation function S* oc (oj) using the Dynamical DMRG (DDMRG) 
[35]. 

2. DMRG results 

The DMRG technique is best suited for open one-dimensional chain sys- 
tems but can be applied to systems with an arbitrary geometry. The icosido- 
decahedron can be viewed as a two-dimensional lattice on a sphere (similar 
to the Kagome lattice, see [8]), i.e., with periodic boundary conditions. This 
means that the convergence is much slower than in one-dimensional systems 
[33]. But since DMRG is a variational method, it is clear that the ground 
state - or the lowest energy in a subspace - is the better the lower the cor- 
responding energy is. Also, the truncated weight Aw offers the possibility 
to judge the quality of the results and an extrapolation to zero truncated 
weight (or m — > oo) might be possible. 

For the investigations throughout this article we employ the Heisenberg 
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Hamiltonian 

¥ = J Yl Ei'£i + 9HBB Si (!) 

i<j i 

with antiferromagnetic isotropic nearest neighbor exchange J only. 

2.1. Numbering of the spins 

When DMRG is applied to spin systems that are not one-dimensional, 
the usual way is to map the system on a one-dimensional chain with long- 
range interactions, i.e., to number the spins of the lattice [36]. However, if 
not very simple systems such as, e.g., ladders are investigated, it is not clear, 
which numbering is best suited. Such a problem also occurs, when DMRG 
is applied in the context of quantum chemistry, where models similar to the 
Hubbard model with long-range interactions appear and the ordering, i.e., 
the numbering of the orbitals is relevant [37, 38, 39, 40, 41]. Since long-range 
interactions diminish the accuracy of DMRG (cf. Ref. [42]) it is clear that a 
good ordering needs to minimize such long-range interactions. 

We have tested several numberings for the icosidodecahedron. The result- 
ing coupling matrices Jij are shown in Fig. 2. The numbering used by Exler 
and Schnack in an earlier investigation [33] (see top left of Fig. 2) gives a very 
regular "interaction pattern" with rather-short-ranged interactions, but the 
"periodic boundaries", i.e., interactions between the first and the last spins, 
are clearly not optimal for the DMRG algorithm with two center sites. As 
proposed in Ref. [37], we have used a variant of the reverse Cuthill-McKee 
algorithm [43, 44], the RCMD algorithm, which aims to number the vertices 
such that the bandwidth of the matrix is minimized. We have also used the 
Sloan algorithm [45] which minimizes the "envelope size", i.e., the sum of 
the "row bandwidths" . (The bandwidth is the maximum of the row band- 
widths.) We have used these algorithms as implemented in Mathematica 
[46]. The figure also shows an unoptimized numbering for comparison. 

The results of DMRG calculations (using the ALPS DMRG code [47]) 
for the different spin numberings are shown in Fig. 3. We have calculated 
the ground state energy of the s — 1/2 icosidodecahedron with an increasing 
number of kept density matrix eigenstates (m) so that the convergence can 
be investigated and a comparison with the exact ground state energy (see 
Ref. [8]) is possible. One can see that the different optimized numberings 
(Exler/Schnack, RCMD, and Sloan) give almost identical results whereas 
the unoptimized numbering gives much worse results. These results show 
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Figure 2: Depiction of the coupling matrix (Jij) for four different numbcrings of the 
vertices of the icosidodecahedron. Nonzero entries are denoted by the orange squares. Top 
left: numbering according to Exler and Schnack (see Ref. [33]); top right: result of the 
RCMD algorithm; bottom left: result of the Sloan algorithm; bottom right: unoptimized 
numbering. 

that a "good" numbering of the spins is absolutely essential if the DMRG 
method is applied to a spin system with a complicated structure. For the 
following results we have always used the numbering as proposed by Exler 
and Schnack. 
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Figure 3: DMRG results for different numberings (see Fig. 2) of the spins sitting on the 
vertices of the icosidodecahedron. The plot shows the error in the ground state energy 
as obtained by DMRG and as a function of the DMRG sweep. The numbers above the 
symbols denote the number of kept density matrix eigenstates for the sweep. 

2.2. Lowest energy eigenvalues and magnetization curves 

As a next step we have calculated the lowest energies in the M subspaces 
for the icosidodecahedron with s > 1/2 using DMRG. The results for the 
s — 1/2 system already showed that DMRG is able to produce very accurate 
results for this system with relative errors smaller than 10~ 3 . 

Fig. 4 shows the lowest energy eigenvalues in the subspaces of total mag- 
netic quantum number M for the icosidodecahedron with s = 1 and s = 3/2 
as obtained by DMRG and - for the large-M subspaces (M > 18 for s = 1 
and M > 33 for s = 3/2) - Lanczos calculations. We have used up to 
m = 2500 density matrix eigenstates for the s = 1 case and up to m = 2000 
for the s = 3/2 case. The largest truncated weight within a sweep is of the 
order of 7 • 10~ 4 for the M = subspace of the s = 1 icosidodecahedron 
and of the order of 4 • 10 -4 for the s = 3/2 case. That the truncated weight 
for the s = 1 icosidodecahedron is larger than for s = 3/2 although more 
states have been used for s = 1 indicates that it cannot be reliably used for 
a quantitative estimate of the error. The reason for this behavior might be 
that the results are not yet fully converged for the value of m that we have 
used, although we have carried out up to 60 sweeps for the calculations. 

The rotational band model predicts a behavior of the form E m i n (M) = 
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Figure 4: Lowest energy eigenvalues in the subspaces of total magnetic quantum number 
M as obtained by DMRG calculations. The eigenvalues for the smallest subspaces, i.e., 
for large M, were calculated using the Lanczos algorithm (calculations performed by J. 
Schnack, private communication). For the DMRG calculations, the ALPS DMRG code 
was employed [47] . For the s = 1 system we have kept m = 2500 density matrix eigenstates 
in all DMRG calculations. For the s = 3/2 system we have kept m = 2000 states for the 
subspaces up to M = 4, m = 1500 states for the subspaces 5 < M < 23, and 1000 states 
for the subspaces M > 23. 

aM(M + 1) + b, i.e., a parabolic dependence [12]. The insets of Fig. 4 
show that this is a good approximation for the energy eigenvalues of the 
full Heisenberg model. The simple rotational band approximation predicts a 
proportionality constant of a = 0.1. The linear fits as shown in the insets give 
the results a = 0.111 for s = 1 and a = 0.108 for s = 3/2, very close to the 
simple rotational band approximation. However, if one uses these (DMRG) 
data to calculate the zero-temperature magnetization curve, it becomes clear 
that there are some crucial deviations from the ideal parabolic dependence. If 
there was an ideal parabolic dependence, the resulting magnetization curve 
would consist of steps with constant widths. Fig. 5 shows the resulting 
zero-temperature magnetization curves as calculated using the DMRG data. 
Again, the exact diagonalization data for s = 1/2 is taken from Ref. [8]. 

One can see, that the magnetization curves do not consist of steps with 
constant widths. There are some anomalies as expected for frustrated sys- 
tems. The plateaus at M./ M. S3it = 1/3 are clearly visible. The magnetization 
jumps due to the independent magnons [26] are also visible. Since the jump 
has a height of AM = 3, it is clear that this effect vanishes for s — > oo 
in the plot because the magnetization is normalized by the saturation mag- 
netization A^sat = SOgfiBS. For s — >■ oo, the classical result, i.e., a strictly 
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Figure 5: Zero-temperature magnetization curves of the icosidodecahedron as obtained 
by DMRG (s = 1, s = 3/2) calculations (cf. Fig. 4). The data for s = 1/2 is taken 
from Ref. [8]. The plot also shows the classical result [19]. The data is normalized to the 
saturation field and magnetization. 



linear magnetization curve [19], will be reached. On the contrary, the plateau 
seems to be very stable when increasing s. A similar behavior has already 
been found for the cuboctahedron, compare Ref. [8]. 

For the cases s = 2 and s = 5/2 we have calculated the lowest energy 
eigenvalues only in some M subspaces, including those subspaces that are 
relevant for the calculation of the plateau width. We have kept m = 2000 
density matrix eigenstates for these calculations. The plateau is clearly vis- 
ible, even for s = 5/2. However, it is not clear if these plateaus can be 
measured in experiments, since for s = 1/2 the plateau disappears rather 
quickly with increasing temperature [21]. Fig. 7 shows the plateau widths 
for different s as obtained by the DMRG calculations (cf. Fig. 6) and the 
exact diagonalization for s = 1/2 [8]. The left part of Fig. 7 shows an ex- 
trapolation to s — > oo using the data for s > 1/2. The extrapolated value for 
the plateau width is 0.013 and not exactly zero but rather close the expected 
value. It has to be emphasized here that the plateau widths shown in Fig. 7 
are approximate values for s > 1/2. The accuracy of the results is analyzed 
in the next subsection. 
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Figure 6: Zero-temperature magnetization curves of the icosidodecahedron for s = 
3/2,2,5/2 as obtained by DMRG calculations (cf. Fig. 4). The plot also shows the 
classical result [19]. The data is normalized to the saturation field and magnetization. 
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Figure 7: Plateau widths for different s quantum numbers as obtained by DMRG S > 1/2 
and exact diagonalization s = 1/2 calculations. The exact diagonalization values are taken 
from Ref. [8]. This part of the plot also includes an extrapolation to s — > oo using the 
data for s > 1/2. The extrapolated value for the plateau width is 0.013. 
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2.3. Extrapolation to m — > oo and error estimates 

A sensible extrapolation in m (the number of kept density matrix eigen- 
states) or Aw (truncated weight) is only possible when the results are fully 
converged for the current m value. Fig. 8 shows the ground state energy and 
the lowest energy in the M = 1 subspace of the s = 1 icosidodecahedron as 
calculated using DMRG for different m values. 
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Figure 8: Left: DMRG results for the lowest energy eigenvalues of the M = and the 
M = 1 subspace. m denotes the number of kept density matrix eigenstates. The lines are 
exponential fits to the data. Right: estimated absolute error as a function of m. 

An extrapolation with an exponential ansatz yields E^* a (M = 0) = 
-41.97(2) J and £g£*(M = 1) = -41.76(1) J. 1 In Ref. [33], the results were 
extrapolated using an ansatz of the form a + b/m. However, this ansatz did 
not work well for our data which is based on much larger m values. The 
right part of Fig. 8 shows the estimated absolute error, i.e., the difference 
between the extrapolated energy and the energy value obtained for finite m, 
as a function of m. For s > 1 such an extrapolation was not possible since the 
production of enough converged energies for different m values is virtually 
impossible. 

Fig. 9 shows the convergence of the width of the .M sa t/3 plateau as a 
function of m. One can see that the plateau width decreases with increasing 
m so that DMRG seems to overestimate the width of plateau. We find similar 
effects for s > 1, but since the calculations are extremely time-consuming, it 
was not possible to obtain enough numerical data for a more systematic study 



1 The errors given in parenthesis are the standard errors of the fitting procedure. 
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Figure 9: The plateau width of the s = 1 icosidodccahcdron as a function of the number 
of density matrix eigenstates that were kept in the DMRG calculations. 



of this effect. Looking more carefully at the truncated weights, we find that 
the DMRG calculations of the lowest energy in the M = 10s(= -M sat /(3g//B)) 
subspace result in smaller values of the truncated weight than the calculations 
in the adjacent subspaces with the same m. However, the order of magnitude 
of the truncated weight is still the same. This indicates that DMRG leads 
to more accurate results exactly at one-third of the saturation magnetization 
which has the consequence that the method seems to systematically overes- 
timate the plateau width if one works with the same m value for adjacent 
M subspaces. However, we estimate that the relative errors in the plateau 
widths are not larger than approximately 10%. 

For many one-dimensional systems DMRG can be considered as a numer- 
ically exact method. This is clearly not the case for the icosidodecahedron. 
Nevertheless, the width of magnetization steps would be accurate if the er- 
rors of the energy eigenvalues are approximately the same for adjacent M 
subspaces. But since the subspace dimensions become smaller for increas- 
ing M it is clear that calculations with fixed m are more accurate for the 
large-M subspaces. However, we suppose that the qualitative features of the 
magnetization curves as predicted by our DMRG calculations are not altered 
by these considerations since the order of magnitude of the errors is about 
the same for energy eigenvalues in adjacent magnetization subspaces. 
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2.4- Comparison with CPS and previous DMRG results 

In this subsection we focus on the s = 5/2 case and the comparison to 
previous results on this system. Fig. 10 shows the ground state energy as a 
function of the DMRG sweep and as a function of the kept density matrix 
eigenstates. One can see that even for 2000 kept states and more than 30 
sweeps the result is not yet converged. A much larger number of states is 
needed to get convergence. Also, an extrapolation to m — > oo is not reliably 
possible because for that many more sweeps would have to be performed. 




sweep # 



Figure 10: Ground state energy of the s = 5/2 icosidodecahedron as a function of the 
DMRG sweep. The numbers show the retained density matrix eigenstates for the current 
sweep. 

For m = 2000 we obtain the value Eq ~ —216.5 J. This value can 
be compared with previous results. The DMRG result of Exler and Schnack 
for the ground state energy (with m = 120) is approximately —211.1 J [33], 
a value that is much higher and thus much more imprecise than our result. 
The very recent result of Neuscamman and Chan using correlator product 
states in combination with variational Monte Carlo is —216.3 J [34]. This 
value is also higher than our DMRG result. Also, the comparison of the 
lowest energies in the M subspaces that are relevant for the calculation of 
the plateau width shows that DMRG is - at least for the CPS ansatz used 
in Ref. [34] - more accurate if enough states are used (see Tab. 1). 
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M 



E DMRG (M)/J 



£ CPS (M)/J 



24 
25 
26 



-158.43 
-153.78 
-147.91 



-154.42 
-149.76 
-144.40 



Table 1: Comparison of the DMRG energies to the CPS energies. For the DMRG calcu- 
lations, m = 2000 states were kept. The CPS data is taken from Ref. [34]. 

2.5. Dynamical correlation function for s = 1/2 

In the following we investigate the dynamical spin correlations of the 
s = 1/2 icosidodecahedron in zero external magnetic field. We focus on the 
local spin autocorrelation function 



found in a non-degenerate, symmetric representation. 2 

We calculate this quantity using two different techniques: i) a contin- 
ued fraction technique based on the Lanczos Exact Diagonalization (ED) 
algorithm [48] and ii) a Dynamical DMRG (DDMRG) approach [35]. 

In Fig. 11 we display the results for S* oc (oj) obtained using both tech- 
niques, with excellent agreement for the same broadening 7] = 0.15J. The 
ED results have been obtained using 500 iterations in the continued fraction 
expansion, while the DDMRG results have been obtained using m = 2000 
states, which results in very time-consuming calculations. In the DDMRG 
f] = 0.15 J is fixed during the calculation, while the ED data can replotted 
using an arbitrary value for 77. For completeness we have also added the pole 
content obtained within the 500 ED iterations. 

Turning to the physics of the spectral function, it displays a rather sharp 
peak at approximately 0.35 J and a very broad shoulder which then falls off at 
approximately 2 J. The rotational band approximation predicts a two-peak 
structure of the dynamical correlation function. This is clearly not the case 
since the "real" spectrum is much broader and shows only one distinct peak. 
However, for such a large s = 1/2 system, we cannot expect the rotational 



2 This is different in some of the ground states realized at finite magnetization, see 




Ref. [8]. 
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Figure 11: The dynamical correlation function Sf oc (uj) as obtained by ED and DDMRG 
calculations. In the ED we have used 500 continued fraction iterations whereas in the 
DDMRG we have used m = 2000 and a broadening 77 = 0.15J. The agreement between 
both methods is perfect. The vertical lines show the energies where the simple rotational 
band model predicts the peaks. 

band model to be a good approximation. This approximation is most accu- 
rate for small systems with large spin quantum numbers s, compare Ref. [49]. 
On the contrary, the strong frustration present in the icosidodecahedron re- 
sults in rather dense spectra of low-lying singlets (even below the lowest 
triplet), low-lying triplets (even below the lowest quintet) and so on [25, 8] 
of which the smeared-out dynamical correlation function is a consequence. 
This behavior is in accord with the experimental INS cross sections obtained 
for the Keplerate Mo72Fe3o with single spin quantum number s = 5/2 [16]. 
It would be very appealing to evaluate the dynamical correlation function 
for this compound, unfortunately this is virtually impossible at the moment. 
The main problem is that rather accurate lowest (S = 0) states need to be 
calculated for the DDMRG method. For the s — 1 icosidodecahedron, the 
estimated error of the ground state energy is already of the order of 0.1J 
using m = 2000. This inaccuracy would be too large for a broadening of 
77 = 0.15 J and either a larger broadening would have to be chosen or more 
states would be needed. However, keeping more states strongly increases 
the calculation time and a much larger broadening would probably blur the 
spectrum too much. Here, a fully optimized and parallelized DDMRG code 
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would be needed that in addition had to run on a supercomputer. 
3. Summary 

In this article we report up to date quantum mechanical DMRG calcu- 
lations for a class of Keplerate magnetic molecules with s = 1/2, .. . ,5/2. 
We demonstrate that it is possible to obtain lowest energy eigenvalues in or- 
thogonal subspaces and magnetization curves for T = with unprecedented 
accuracy. In addition, using DDMRG as well as Lanczos techniques, we ob- 
tained the dynamical correlation function for Keplerates with an intrinsic 
spin of s = 1/2. 
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